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J..   INTRODUCTION 

In  1962  Gene  Golub  questioned  why  I  was  starting  to  work  on  the 
numerical  ODE  problem.   "Haven't  all  the  important  problems  been  solved?" 
he  asked.   Indeed,  it  seemed  at  the  time  that  they  had.   Dahlquist's 
important  paper  [8]  had  appeared  6  years  earlier  and  solved  what  seemed  to 
be  all  of  the  important  and  interesting  theoretical  questions;  most 
numerical  analysts  were  unaware  of  any  major  classes  of  ordinary 
differential  equations  that  could  not  be  solved  by  computer  at  speeds  that 
seemed  to  be  reasonable  for  the  difficulty  of  the  problem;  and  the 
available  computer  codes  seemed  to  perform  fairly  well.   At  the  most  it 
appeared  that  there  were  a  few  "clean-up"  questions  to  take  care  of. 
However,  at  that  time,  the  problem  of  stiffness  was  already  known  to  some 
people.   It  appears  to  have  first  appeared  in  1952  (Curtiss  and 
Hirschfelder  [7])  and  Dahlquist  was  about  to  publish  his  classic  1963  paper 
[9]  on  it.   Ten  years  later,  reasonably  robust  software  began  to  appear, 
but  even  today  there  are  still  several  important  theoretical  questions 
waiting  to  be  answered  that  may  have  significant  practical  implications  for 
stiff  problems. 
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Apart  from  these  questions,  is  there  anything  else  to  do  in  the  field? 
This  paper  examines  the  steps  in  the  development  of  a  typical  problem  from 
the  time  it  is  first  identified  by  the  "user"  (that  is,  by  the 
scientist/engineer  who  runs  into  an  intractable  computational  problem)  to 
the  time  that  robust  software  is  developed  for  classes  of  problems.   This 
process  may  initially  take  20  years;  during  the  early  years  very  few  will 
be  aware  of  the  problem,  and,  even  at  the  end  after  robust  software  has 
appeared,  development  may  continue.   After  examination  of  the  development 
process,  the  paper  looks  at  a  number  of  problems  that  are  in  various  stages 
of  development  and  reviews  recent  results  in  these  areas .   We  will  see  that 
there  are  a  number  of  problems  of  great  interest  at  all  stages  of 
development  and  that  not  only  have  there  been  a  number  of  important  new 
results  recently,  but  that  there  are  many  open  problems. 

_1  ._1.   The  Development  Staircase 

The  progress  from  problem  identification  to  robust  software  can  be 
broken  into  six  levels.   These  are: 

1  Identification  of  problem  by  a  user. 

2  Development  of  a  method  and  preliminary  codes. 

3  Mathematical  analysis  of  the  problem  and  method. 

4  Preparation  of  initial  production  codes. 

5  Analysis  of  software  for  classes  of  problems. 

6  General  purpose,  robust  software  development. 
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In  the  case  of  stiff  ODEs,  these  developments  took  two  decades  from  the 
early  1950's.   Today,  the  solution  of  most  stiff  problems  is  relatively 
routine  unless  additional  problems  (such  as  very  large  systems  or  many 
discontinuities)  are  present.   However,  the  development  has  not  stopped. 
New  methods  and  new  analyses  continue  to  appear,  giving  rise  to  the 
"staircase"  picture  shown  in  Figure  1.   Each  level  of  the  staircase 
continues  to  spawn  new  developments  at  higher  levels. 

The  stiff  problem  has  already  reached  the  top  of  this  staircase,  so, 
once  again,  we  might  ask  "Is  there  anything  left  to  do?"  besides  a  few 
clean-up  questions.   This  paper  examines  a  number  of  problems,  some  of 
which  are  at  or  near  the  top  of  the  staircase,  but  others  of  which  are  only 
on  the  first  step.   Variable-order,  variable-step  methods,  for  example  (or, 
more  precisely,  problems  for  which  variable-order,  variable-step  methods 
are  essential  for  efficient  integration)  are  also  on  the  top  level,  but 
problems  which  require  low  accuracy  methods  are  at  the  first  step:  they 
have  been  identified  but  little  has  been  done  for  their  solution.   The 
problem  areas  we  will  examine  are : 

Stiff  problems 

Variable-order,  variable-step  methods 

Problems  with  eigenvalues  near  the  imaginary  axis 

Oscillating  problems 

Very  large  stiff  problems 
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Problems  with  frequent  discontinuities 

Multirate  methods 

Low  accuracy  methods 

The  approximate  status  of  these  problems  is  shown  in  Figure  2.   This  figure 
does  not  include  any  of  the  problems  in  ODEs  arising  from  PDEs  such  as 
variable-sized  systems  that  occur  when  the  method  of  lines  is  used  with  a 
varying  space  discretization,  and  the  special  structure  of  the  Jacobian 
that  arises  from  the  method  of  lines  which  may  make  very  large  problems 
more  tractable  than  in  the  general  case.   If  we  had  a  crystal  ball,  we 
could  extend  Figure  2  to  see  what  problems  we  will  be  facing  in  the 
future — but  that  would  automatically  put  those  problems  on  step  1  by  virtue 
of  their  having  been  identified!   However,  the  fact  that  there  are  problems 
at  all  levels  suggests  that  there  will  no  shortage  of  problems  in  the 
future. 

_2.   Problems  showing  recent  progress 

This  section  examines  some  results  of  the  last  few  years  that  have 
contributed  to  our  understanding  of  some  of  the  problems  and  methods  for 
solving  them.   The  results  mentioned  are  frequently  applicable  to  several 
of  the  problems;  in  the  space  available  we  will  discuss  their  application 
only  to  one  problem. 

2«_1.   Stiff  problems 

The  major  outstanding  problem  of  the  last  few  years  bas  been  the  error 
and  stability  analysis  of  stiff,  non-linear  systems.   Burrage  and  Butcher 
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[4]  and  Burrage  [3]  have  examined  the  stability  of  non-linear,  non- 
autonomous  systems  for  Runge-Kutta  methods.   They  have  derived  a  sufficient 
algebraic  condition  (i.e.,  a  condition  on  the  coefficients  of  the  method) 
for  stability  of  the  numerical  solution  to  problems  which  satisfy 

(2.1)  Re<f(y,t)  -  f(z,t),y  -  z>  <  0 

Note  that  this  condition  implies  the  stability  of  the  differential  system. 
Stability  of  the  method  is  defined  to  be  the  condition 

(2.2)  Hy  -  z  ||  <  ||y    -  z    || 

n    n        n—  1    n-l 

This  condition  implies  that  the  difference  between  numerical  solutions 

reduces  at  each  step.   A  similar  approach  is  now  being  used  for  multistep 

method  analysis  by  Dahlquist  [10],  [11],  and  Nevalinna  and  Liniger  [19]  who 

use  eq.  (2.2)  with  y  and  z  interpreted  as  vectors  of  all  past  values  used 

in  computing  the  forward  solution.   This  is  a  much  stronger  condition  than 

absolute  stability  which  only  requires  that  | |y  -  z   |  be  bounded  as  n 

n    n 

goes  to  infinity.   The  former  condition  is  called  contractivity.   The 
advantage  of  working  with  contractivity  is  that  it  is  possible  to  examine 
each  step  independently.   The  difficulty  with  the  analysis  of  stability  of 
multistep  methods  is  that  adjacent  steps  interact,  so  that  non-linearities 
in  the  problem  or  variable  stepsizes  give  rise  to  products  of  changing 
operators  for  which  there  may  not  be  a  single  fixed  norm  in  which  they  are 
all  bounded  by  one  although  the  product  may  be  bounded .  When  we  have 
contractivity  in  some  fixed  norm  for  the  range  of  problems  and  methods 
used,  it  automatically  becomes  a  global  property. 

The  analysis  has  been  greatly  aided  by  the  introduction  of  one-leg 
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methods  by  Dahlquist  (see  [10]).   These  are  intimately  related  to  multistep 
methods,  taking  the  form 

k  k 

(2.3)      I   ay    +  hf(  I    By   )  -  0 
1-0  x   n'i      i=0  x   n_i 

In  fact,  for  fixed-step  methods,  these  methods  are  equivalent  to  multistep 

methods  with  the  same  coefficients.   However,  they  differ  for  variable-step 

methods.   The  one-leg  methods  are  easier  to  analyse  and  probably  have 

better  stability  properties.   The  ease  of  analysis  arises  from  the  fact 

that  each  step  is  governed  by  a  difference  formula  in  which  only  one  value 

of  the  function  f  is  computed.   Contractivity  can  be  studied  by  local 

linearization  of  the  problem  and  there  is  only  one  Jacobian  to  deal  with  in 

a  linearization  of  a  one-leg  method.   For  example,  if  we  apply  the 

trapezoidal  rule  to  stepsize  h  ,  .  /0  =  t  ,  .  -  t  ,  namely, 

n+l/z    n+1    n 

y  .i  =  y  +h  ,,  ,„(f  ,.  +  f  )/2,  we  find  that  contractivity  requires 
n+1    nn+l/2n+ln 

II"  -k+l/2Jn+l1"'tI+it,n+l/2Jn]ll  '    l 

where  J  and  J  .,  are  the  Jacobians  at  T  and  t  ,,  respectively.   Since  J 
n      n+1  n      n+1 

may  change  over  the  interval,  this  bound  may  be  very  difficult  to 

establish;  indeed  it  may  not  be  true.   On  the  other  hand,  the  equivalent 

one-leg  method,  namely  y  ,,  =  y  +  h  ,,  y„f((y  .,  +  y  x  /0\  leads  to  a 

n+1    n    n+1/2    n+1    n)/2) 

requirement  that 

i|[I--k+i/2jrl(I+i\,+i/2J!"  *  > 

where  J  is  evaluated  at  t  . .  ,_  in  both  instances  in  the  last  equation. 

n+1/2 

Dahlquist  [10]  has  shown  that  if  the  problem  satisfies  eq.  (2.1),  a 
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fixed  step,  constant  coefficient,  A-stahle  multistep  method  is  such  that 
tbere  exists  a  norm  in  which  the  method  is  contractive  and  conversely, 
''ence,  A-stahility  is  an  attractive  property  for  multistep  methods. 
Unfortunately,  A  stability  limits  the  order  to  two  so  there  is  considerable 
interest  in  generating  similar  results  for  less  stringent  stability 
conditions  such  as  A(a)  and  stiff  stability.   Somewhat  surprisingly,  the 
algebraic  stability  condition  of  Burrage  and  Butcher  [4]  is  stronger  than 
A-stabilitv,  meaning  that  an  A-stable  PK  method  may  not  be  contractive  for 
a  problem  satisfying  eq.  (2.1).   Since  contractivity  and  stability  of  one- 
step  methods  are  essentially  identical  (a  method  that  is  not  contractive  in 
any  steps  is  hound  to  he  unstable),  an  A-stable  RK  method  may  not  be 
suitable  for  a  non-linear  stiff  equation. 

_?._?.   Problems  with  imaginary  eigenvalues 

If  the  system  Jacohian  has  eigenvalues  near  the  imaginary  axis  and  the 
corresponding  components  are  known  to  be  absent  from  the  solution,  one  is 
interested  in  a  method  which  damps  such  components  numerically, 
independently  of  the  stepsize  used.   Since  absolute  stability  is  desired 
for  eigenvalues  with  more  negative  real  parts,  one  is  forced  to  ask  for  A- 
stability  (unless  a  special  nethod  that  does  exponential  fitting  is  used). 
'Hie  limitation  on  order  of  multistep  methods  that  are  A-stahle  forces 
consideration  of  other  methods.   Two  sets  of  recent  results  are 
particularly  significant  here.   A  charminp  paper  by  banner,  wairer,  and 
Norsett  f?5]  answers  many  open  questions  about  the  maximum  order  of  a 
variety  of  A-stable  methods  and  gives  simple  proofs  for  many  Vnown  results. 
In  narticular,  it  is  shown  that  the  maximum  order  of  an  A-stable  method 
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cannot  exceed  2q,  where  q  is  the  number  of  derivatives  or  stages  used. 
This  means  that  the  choice  for  higher-order  A-stable  methods  is  restricted 
to  highly  implicit  methods  such  as  implicit  Runge-Kutta  methods  and  multi- 
derivative  methods.   The  latter  are  not  very  attractive  because  the 
Jacobian  must  be  computed  exactly.   (They  may  be  suitable  for  problems  in 
which  the  Jacobian  is  known  explicitly  or  is  fixed.)   Implicit  Runge-Kutta 
methods  are  potential  candidates,  but  the  size  of  the  non-linear  problem 
that  has  to  be  solved  is  a  problem.   For  systems  of  s  equations  using  a  q- 
stage  method,  a  simultaneous  system  of  sq  equations  must  be  solved.   Work 
by  Bickart  [2]  and  Butcher  [6]  shows  how  to  implement  implicit  Runge-Kutta 
methods  in  a  relatively  efficient  manner  but  a  more  promising  approach  is 
to  use  semi-implicit  Runge-Kutta  methods  with  equal  diagonal  elements  so 
that  only  one  LU  decomposition  is  needed  along  with  q  back  substitutions 
for  a  q-stage  method.   Unfortunately,  these  q-stage  methods  are  limited  to 
order  q+1.   Another  approach  being  explored  by  Butcher  [5]  is  based  on 
iterated  defect  correction,  and  also  requires  the  solution  of  a  set  of  s 
non-linear  equations.   Yet  another  scheme  is  the  blended  multistep  method 
of  Skeel  and  Kong  [23].   This  is  a  linear  combination  of  two  multistep 
methods,  where  the  combination  depends  on  the  Jacobian.   The  accuracy  is 
derived  from  the  methods  and  is  not  dependent  on  the  accuracy  of  the 
Jacobian.   On  the  other  hand,  if  the  stability  for  linear  problems  is 
analyzed  using  the  true  Jacobian,  the  combined  method  can  have  the 
stability  properties  of  a  multiderivative  method.   The  multiderivative  and 
equivalent  blended  methods  require  the  solution  of  a  linear  equation 
involving  a  polynomial  in  the  Jacobian.   If  this  polynomial  is  chosen  to 
have  equal  roots,  an  LU  factorization  of  a  linear  polynomial  can  be 
performed,  followed  by  multiple  back  substitutions.   This  makes  maximum  use 
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of  the  sparsity  by  avoiding  powers  of  the  Jacobian. 

2.3.      Oscillating  Problems 

We  distinguish  between  problems  with  almost  purely  imaginary 
eigenvalues  corresponding  to  components  not  present  in  the  solution  and 
problems  with  oscillating  solutions.   In  the  former,  the  problem  is  to 
avoid  introducing  the  oscillating  components.   Since  difficulties  arise  in 
the  former  case  only  when  the  imaginary  parts  are  very  large,  these  have 
been  viewed  in  the  same  category  as  stiff  problems,  and  this  seems  to  be  a 
reasonable  viewpoint.   (See  Miranker  and  Wahba  [18].)  If,  however,  the  true 
solution  exhibits  very  high  frequency  oscillations,  the  problem  can  be  very 
different.   Only  if  the  effect  of  the  oscillations  is  linear  is  it 
reasonable  to  damp  them  and  seek  the  average  solution.   In  this  restricted 
case,  it  is  still  reasonable  to  view  the  problem  as  a  stiff  one  because  the 
objective  of  the  numerical  method  is  the  same,  _to  ignore  components  due  to 
large  eigenvalues .   However,  if  the  effect  of  the  oscillation  on  the  rest 
of  the  system  is  non-linear,  it  is  not  possible  to  damp  the  oscillation 
without  error.   In  this  more  general  case,  it  is  not  reasonable  to  view  it 
as  a  stiff  problem;  rather  we  must  ask  what  it  is  we  want  to  determine 
about  the  solution.   Some  recent  work  by  Petzold  and  Gear  [22]  and  Petzold 
[21]  addresses  this  issue. 

Suppose  we  want  to  know  the  waveform  after  a  very  large  number  of 
cycles,  or  the  envelope  of  the  solution  (see  Figure  3).   Unfortunately,  the 
waveform  is  not  defined  for  an  initial  value  problem  because  the  initial 
value  defines  a  unique  solution.   However,  the  "envelope"  in  Figure  3  is 
apparent  to  any  engineer!   To  make  the  problem  mathematically  (and 
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numerically)  tractable,  we  define  a  quasi  envelope.   This  is  a  smooth  curve 
which  can  be  "integrated"  using  stepsizes  which  may  be  very  much  larger 
than  the  period  of  the  oscillation  T.   The  definition  of  the  quasi  envelope 
is  illustrated  in  Figure  A.   If  the  period  T  is  specified  for  an  autonomous 
differential  equation,  the  quasi-envelope  z(t)  is  defined  in  terms  of  an 
arbitrary  initial  segment  for  t  in  [0,T)  by  integration  over  single 
periods.   Thus,  the  value  of  z  at  s+t  is  defined  in  terms  of  z(s)  by  the 
formula 

z(s+T)  =  z(s)  +  (v(s+T)  -  v(s)) 

where  v  is  the  solution  of  the  differential  equation  with  initial  values 
(s,z(s)).   Clearly,  if  z(0)  =  y(0),  z  agrees  with  the  solution  of  the 
initial  problem  at  all  integral  multiples  of  T.  We  will  write  the  above 
equation  in  the  form 

(2. A)     z(t+T)  =  z(t)  +  Tg(t,z(t)) 

where 

(2.5)     g(t,z)  =  (v(t+T)  -  v(t))/T 

Equation  (2. A)  looks  like  the  Euler  method  applied  to  the  differential 
equation  z'  =  g(t,z),  and  this  idea  is  exploited  in  developing 
interpolation  methods.   After  some  manipulation  we  can  devise 
approximations  that  look  like  integration  methods  in  which  the  coefficients 
are  functions  of  T/H  where  H  is  the  stepsize  used.   As  this  ratio  tends  to 
zero,  they  tend  to  standard  integration  techniques.   (Note  that  equation 
(2. A)  tends  to  a  theN differential  equation  z'  =  g(t,z)  as  T  tends  to  zero.) 
For  example,  the  multistep  methods  take  the  form 
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k 

(2.6)      I    (a  z    +   H0  g    )  =  0 
j_n    n— 1     1  n— i 

If  this  is  exact  for  polynomials  z  of  degree  k  in  t,  the  coefficients  a 
and  3  are  polynomials  of  degree  at  most  k  in  T/H.   Note  that  each  step  of 
this  "multistep  method"  requires  the  evaluation  of  g(z),  and  this  requires 
the  integration  of  the  original  differential  equation  through  one  cycle. 

If  the  period  T  is  unknown,  or  if  it  varies  slowly,  it  can  be 
estimated  by  minimizing 

t+T 
/  ll(y(s+T)  -  y(s)irds 
t 

with  respect  to  T.   In  fact,  this  serves  to  define  T(t)  since  it  is  not 

otherwise  defined.   The  previous  formulas  were  stated  for  a  fixed  period  T. 

The  same  formulas  can  be  used  for  a  variable  period  by  using  a  change  of 

independent  variable  to  one  in  which  the  period  is  constant .   When  the 

equations  are  expressed  in  the  new  independent  variable,  we  get  a 

relationship  similar  to  eq.  (2.4)  but  with  a  fixed  period.   In  addition,  we 

must  add  an  equation  for  the  real  time.   It  takes  the  form 

t(x  +  u)  =  t(x)  +  T(t) 

where  jj  is  the  fixed  period  in  the  new  independent  variable  t. 

The  general  form  of  the  method  consists  of  a  pair  of  integrators,  the 
inner  integrator  called  from  a  program  which  determines  g  and  estimates  T, 
plus  the  outer  integrator  which  computes  z.   Each  step  of  the  outer 
integrator  requires  one  or  more  evaluations  of  g  in  a  typical  predictor- 
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corrector  operation.   In  fact,  the  outer  integration  can  be  done  by  minor 
modifications  of  a  conventional  integrator,  and  the  step  and  order 
selection  algorithms  can  be  carried  over.   Each  time  the  the  function  g  is 
evaluated  for  the  outer  integration,  an  integration  over  a  full  cycle  by  a 
conventional  integrator  is  performed.   If  the  period  can  vary,  it  is 
necessary  to  integrate  over  two  full  cycles  to  compute  the  term  which  has 
to  be  minimized. 

Non-autonomous  problems  in  which  the  t-dependency  of  the  differential 
equation  contains  frequencies  which  determine  the  period  can  be  handled  by 
the  same  techniques,  but  the  large  stepsizes  H  must  be  synchronized  with 
the  period.   If  the  t-dependencies  are  very  slow  compared  to  the  period, 
this  synchronization  is  not  important. 

The  method  is  not  efficient  if  the  problem  also  happens  to  be  stiff 
(meaning  that  there  are  rapidly  decaying  components  present  also)  because 
it  is  then  necessary  to  use  an  implicit  outer  integration  scheme  which 
means  that  a  Jacobian  must  also  be  calculated.   In  this  case,  the  Jacobian 
is  a  state  transition  matrix  over  one  period,  and  this  matrix  is  never 
sparse  even  if  the  system  is  sparse.   The  method  is  also  not  effective  if 
the  solution  is  composed  of  several  almost  periodic  functions  with  nearby 
periods.   These  types  of  problems  remain  almost  intractable. 

2i'J±.      Discontinuous  problems 

After  each  discontinuity,  earlier  information  cannot  be  used  to 
extrapolate  forward  so  either  a  single-step  method  must  be  used  or  a 
restart  of  a  multistep  method  is  necessary.   Most  current  automatic  codes 
for  multistep  methods  use  as  a  first  step  a  first-  or  second-order  one-step 
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method  that  is  also  one  of  their  multistep  methods,  and  then  use  the 
variable-order  mechanism  to  build  up  the  order.   This  is  moderately 
efficient,  and,  if  it  occurs  infrequently  is  an  adequate  technique. 
However,  frequent  restarting  is  costly.   Therefore,  we  have  developed  a 
Runge-Kutta  starter  for  multistep  methods  that  uses  a  single  Runge-Kutta 
like  step  to  generate  all  the  information  needed  to  switch  to  a  moderately 
high  order  multistep  method  in  the  second  step  [16].   The  method  consists 
of  a  single  step  of  the  form 


k   =  hf(y   +  I    3k)  ,    1-1,  ...,  q, 


ym  =  yn  +  l   Y.k   ,    m  =  1,  ...,  s. 

This  is  a  q-stage  R-K  method  that  attempts  to  generate  s  new  values  in  one 
step.   These  s  values  can  be  used  to  start  a  multistep  method  of  order  s+1 
or  s+2.   It  has  been  shown  that  for  implicit  RK  starters,  it  is  possible  to 
generate  a  set  of  s  points  of  order  s  using  s  stages,  and  for  the  more 
promising  explicit  methods,  we  can  generate  a  set  of  fourth-order  points 
using  6  stages.   Hence,  a  non-stiff  method  can  be  started  at  fourth  order 
using  six  evaluations.   In  fact,  there  is  a  nine-parameter  family  of  such 
methods.   Preliminary  tests  indicate  that  it  reduces  the  number  of  function 
evaluations  to  get  started  by  about  half. 

A  problem  that  still  remains  to  be  solved  is  that  of  detecting  the 
presence  of  a  discontinuity  efficiently.   If  nothing  is  known  about  its 
nature,  a  code  has  great  difficulty  reducing  the  step  until  it  lands  right 
on  the  discontinuity  with  a  mesh  point. 
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3^,      Problems  presenting  major  difficulties 

This  section  examines  some  of  the  problems  which  appear  to  be  more 
difficult  to  solve.   Perhaps  the  most  important  of  these  is  the  large, 
stiff  problem,  both  because  it  is  a  practical  problem,  and  its  solution 
will  have  an  important  impact  on  the  application  of  ODE  methods  to  PDE 
problems  via  the  method  of  lines.   The  major  remaining  problems  with 
variable-step,  variable-order  codes  arise  from  the  desire  to  put  them  on  a 
firm  theoretical  basis.   In  general,  they  work  very  well  and  further 
advances  are  unlikely  to  improve  their  performance  significantly.   The  last 
two  problems,  multirate  methods  and  low  accuracy  methods,  are  in  their 
infancy.   The  latter  is  very  important  in  real-time  simulation  in  which  the 
integration  is  time-critical,  that  is,  it  is  of  no  interest  unless  one 
second  of  simulated  time  can  be  computed  in  no  more  than  one  second  of  real 
time!   This  problem  may  also  be  important  in  the  method  of  lines  for  PDEs 
where  typical  accuracies  (two  digits)  are  very  low  compared  to  typical 
accuracies  requested  in  ODE  solution  (5  to  15  digits) . 

_3._1.   Large  stiff  problems 

We  can  distinguish  two  categories  of  "largeness"  in  stiff  problems: 
those  in  which  an  LU  decomposition  of  the  Jacobian  is  possible  but  a  time 
consuming  part  of  the  calculation,  and  those  in  which  it  is  impossible  to 
contemplate  factoring  the  Jacobian  because  of  space  considerations.   In  the 

latter  case  it  may  be  impossible  to  store  the  Jacobian  explicitly,  even  if 

3 
it  is  sparse.   The  former  consist  of  order  of  10  equations,  while  the 

latter  may  have  as  high  as  10  equations  or  larger.   The  latter  usually 

arise  from  two-  and  three-dimensional  PDEs  treated  by  the  method  of  lines. 


-  15  - 


If  the  number  of  large  eigenvalues  causing  stiffness  is  small,  it  may 
be  possible  to  treat  the  problem  efficiently  using  ideas  such  as  those  in 
the  Correction  in  the  Dominant  Subspace  (CDS)  method  described  in  Alfeld 
and  Lambert  [1].   The  basic  idea  behind  this  class  of  methods  can  be 
explained  as  follows.   If  a  stiff  method,  such  as  backward  Euler,  is 
applied  to  the  problem 


y'  =  Ay 
we  get  a  difference  equation  of  the  form 


Vl  "   R(hA)yn 


where  R(hA)  =  [I  -  hA]   .   The  effect  of  large  eigenvalues  in  hA  is  to 

remove  the  corresponding  eigenvectors  in  y  ,..   Small  eigenvalues  cause  the 

n+l 

corresponding  eigenvectors  to  be  multiplied  by  an  order  one  approximation 

hA 
to  e   .   If  an  explicit  method  such  as  Euler  is  used,  we  find  that  R(hA)  = 

[I  +  hA].  This  correctly  treats  the  eigenvectors  corresponding  to  the  small 

eigenvalues,  but  amplifies  those  corresponding  to  large  eigenvalues.   The 

CDS  method  consists  of  applying  a  simple  (probably  explicit)  method,  and 

then  damping  the  eigenvectors  corresponding  to  the  large  eigenvalues  (the 

dominant  subspace).   This  requires  a  knowledge  of  the  dominant  subspace 

but,  if  its  dimension  is  small,  does  not  involve  much  storage  or 

calculation.   The  damping  is  applied  by  removing  from  the  solution  most  of 

the  component  in  the  dominant  subspace.   Details  can  be  found  in  the  cited 

paper.   The  method  requires  a  few  inner  products  and  vector  operations  so 

that  the  calculation  is  proportional  to  the  dimension  of  the  system  and  the 

dimension  of  the  dominant  subspace;  much  better  than  the  n  or  n^ 

operations  required  with  a  full  LU  handling  of  the  Jacobian.   The  potential 
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difficulties  with  this  sort  of  approach  are  that  a  time  varying  Jacobian 
means  a  relatively  costly  recalculation  of  the  dominant  supspace,  and  that 
the  vectors  spanning  a  dominant  subspace  are  not  sparse  even  when  the 
matrix  is  very  sparse  (except  in  unusual  circumstances) . 

If  a  conventional  direct  method  is  applied  to  problems  in  the 

intermediate  size  range,  the  linear  equation  subproblem  has  two  important 

3 
components,  the  LU  decomposition  which  may  be  as  slow  as  n  operations,  and 

2 
the  back  substitution  which  may  be  as  slow  as  n  operations.   A  new 

decomposition  must  be  performed  under  a  number  of  conditions.   These  are 

when  the  actual  Jacobian  changes  and  when  the  order  or  stepsize  changes. 

The  latter  occurs  because  the  matrix  which  must  be  decomposed  has  the  form 

[I  -  h3J]  where  J  is  the  system  Jacobian.   Enright  [12]  suggested  one  way 

of  avoiding  a  complete  refactorization  when  only  h  or  3  changes  (the  most 

common  case) .   It  involves  factoring  the  matrix  into  SHS   where  S  is  a 

suitable  similarity  or  unitary  transformation,  and  H  is  an  easy  matrix  to 

work  with  (for  example,  an  upper  Hessenberg  matrix).   Changes  in  h  and  3 

can  then  be  accomodated  by  adding  a  multiple  of  the  identity  to  H.   Enright 

and  Kamel  [13]  have  since  suggested  another  scheme  in  which  S  is  a  unitary 

transformation  (a  product  of  Householder  transformations  for  example),  and 

a  permutation  matrix  is  also  folded  into  S  so  that  all  of  the  large 

elements  of  H  appear  in  the  first  rows.   The  reduction  to  this  form  can  be 

stopped  at  an  intermediate  point  when  H  consists  of  large  elements  in  the 

upper  rows,  zeros  in  the  bottom  left-hand  part,  and  small  elements  in  the 

bottom  right  hand  part.   At  this  point,  the  small  elements  can  be  discarded 

and  the  resulting  matrix  used  as  an  approximation  in  the  iteration.   It 

achieves  an  effect  similar  to  correction  in  the  dominant  subspace  since  it 
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damps  the  dominant  space  by  applying  an  implicit  method  with  a  quasi-Newton 
iteration  to  them.   (Of  course,  the  underlying  method  in  this  scheme  must 
have  appropriate  stability  properties  for  stiff  problems.) 

In  the  case  of  extremely  large  systems,  it  appears  that  the  only  way 
to  handle  the  problem  is  by  iterative  methods  so  that  at  most  the  original 
Jacobian  has  to  be  stored.   (If  the  equations  are  derived  by  the  method  of 
lines,  the  Jacobian  elements  often  can  be  generated  as  needed,  so  even  they 
need  not  be  stored  explicitly.)  Examination  of  the  ideas  used  in  the  CDS 
method  suggests  that  the  important  property  of  any  iterative  method  is  that 
it  must  be  reasonably  accurate  for  components  in  the  dominant  subspace. 
Warming  and  Beam  [26]  explore  approximate  factorizations  by  operator 
splittings  that  are  related  to  ADI  methods.   In  these,  the  equivalent  of 
the  matrix  [I  -  h(3J]  is  approximated  by  the  product  [I  -  h3J,]  [I  -  h$J_] 
where  J  =  J.  +  J„.   The  matrices  J.  and  J„  contain  the  components  due  to 
derivatives  in  the  x  and  y  spatial  directions,  respectively.   For 
rectangular  regions,  the  eigenvalues  of  J  are  the  sums  of  eigenvalues  of  J. 
and  J-.   Further,  the  eigenvectors  corresponding  to  large  eigenvalues  of 
either  J.  or  J„  are  in  the  null  space  of  the  other.   As  long  as  this 
property  is  approximately  true,  the  factorization  is  a  reasonably  accurate 
approximation.   Since  the  factors  correspond  to  one-dimensional  operators, 
they  are  effectively  tridiagonal  for  second-order  differences,  and  suitably 

banded  for  higher-order  differences.   Unfortunately,  terms  involving 

2 
3  /9x3y  do  not  permit  this  factorization,  but  frequently  these  terms  do  not 

contribute  to  the  stiffness  which  arises  from  parabolic  terms. 

Consequently,  they  can  be  ignored  in  the  quasi-Newton  step. 

The  underlying  theme  of  all  of  these  methods  is  similar,  and  can  be 
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related  to  the  use  of  one  or  more  predictor-corrector  schemes  where  the 
corrector  iteration  is  performed  using  quasi-Newton  in  which  the  Jacobian 
is  approximated  by  another  matrix  chosen  to  get  a  "simpler"  factorization. 
The  important  property  of  the  approximation  is  that  it  must  lead  to  fast 
convergence,  and  this  is  determined  by  its  ability  to  represent  the 
Jacobian  reasonably  accurately  in  the  dominant  subspace. 

_3._2.   Variable-step,  variable-order  methods 

An  adequate  theory  for  automatic  methods  is  still  lacking.   A  theory 
is  needed  to  justify  the  error  estimators  used  to  control  the  step  and 
order.   Most  current  estimators  are  based  on  the  asymptotic  analysis  of 
fixed-step,  fixed-order  methods  in  which  a  power  series  in  the  stepsize  h 
is  developed  for  the  error.   Higher-order  terms  are  dropped  and  differences 
of  the  solution  are  used  to  estimates  the  derivatives  that  occur  in  the 
leading  terms.   A  variable-order,  variable-step  method  is  controlled  by  a 
user  parameter,  say,  £,  which  specifies  the  error  tolerance.   If  the  global 
error  and  local  error  were  a  smooth  function  of  £,  one  could  consider  an 
asymptotic  expansion  in  £  instead  of  h,  but  even  if  a  code  attempts  to 
adjust  parameters  so  that  the  local  error  is  exactly  equal  to  the  desired 
multiple  of  £  at  each  step,  the  global  error  would  still  not  be  a 
reasonable  function  of  £  because  the  switch  from  one  order  to  another  is  a 
discrete  process.   In  practice,  codes  do  not  achieve  the  desired  local 
error  exactly,  but  accept  any  step  that  appears  to  give  an  error  no  larger 
than  desired.   Furthermore,  they  tend  to  avoid  too  many  changes  to  stepsize 
and  order  in  the  interest  of  efficiency,  and  this  tends  to  make  the  global 
error  an  erratic  function  of  £.   As  £  is  reduced,  the  global  error 
generated  by  the  code  should  reduce.   Ideally,  the  reduction  should  be 
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proportional,  but  this  seems  to  be  difficult  if  not  impossible  to  achieve. 
See  Stetter  [24].   An  alternative  approach  could  be  based  on  the 
effectiveness  theory  of  Hull  [17].   That  is  based  on  showing  that  a  code 
satisfies  the  error  requirements  for  any  problem  from  a  specified  class  of 
problems.   However,  this  is  not  a  trivial  task,  and  it  is  not  usually 
possible  to  verify  that  the  problem  lies  within  the  given  class.   It  is 
hard  to  see  any  promising  directions  for  this  problem. 

_3._3.   Multirate  methods 

A  multirate  method  is  a  technique  to  handle  each  variable  with  a 
stepsize  appropriate  to  its  behavior.   Thus,  if  the  system  of  equations  is 


dyi 

the  i-th  variable  should  be  integrated  using  a  stepsize  function  h  (t). 
Each  f  must  be  capable  of  being  evaluated  separately  because  the  t  values 
at  which  one  f  must  be  evaluated  may  be  different  from  the  t  values  used 
for  another.   Since  the  f.'s  depend  on  the  y   for  other  values  of  j,  it  is 
necessary  to  interpolate  for  the  values  of  y  at  the  values  of  t  used  as 
mesh  points  for  y.  but  not  for  y  .   The  cost  of  interpolation  is  non- 
negligible  compared  to  the  cost  of  an  integration  step  (in  fact,  it  is 
about  the  cost  of  a  predictor  step) .   Hence,  multirate  methods  can  only  be 
competitive  if  one  of  two  conditions  is  true:  either  the  cost  of  function 
evaluations  must  be  very  high  or  the  dependency  of  one  function  f  on  other 
values  of  y  must  be  sparse.   Then  there  can  be  considerable  savings.   In 
the  former  case  the  savings  can  arise  because  the  total  number  of  f's 
evaluated  will  be  reduced;  in  the  latter  case  the  savings  can  arise  because 
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very  few  interpolations  are  needed.   The  latter  case  occurs  when  the  ODE's 
arise  from  a  method  of  lines  treatment  of  PDE's. 

Even  if  the  necessary  conditions  are  met,  there  are  many  problems  to 
be  solved  before  the  method  is  practical.   The  first  concerns  the  errors 
produced  by  the  interpolation  process.  While  it  is  simple  to  show  that  a 
straightforward  scheme  has  the  appropriate  order  of  convergence  (see  Gear 
[14]),  it  is  necessary  to  derive  estimation  formulas  for  the  errors  created 
by  each  integration  process  so  that  the  stepsizes  can  be  controlled 
automatically.   Perhaps  the  usual  estimates  of  the  derivatives  will 
suffice,  but  there  is  a  danger  that  interpolation  errors  from  other 
components  may  interfere  because  these  errors  will  have  varying  signs  even 
over  the  space  of  a  few  steps  since  they  depend  on  the  relationship  between 
the  mesh  used  in  one  equation  and  that  used  in  another. 

When  we  have  error  estimation  formulas,  we  have  to  determine  how  to 
structure  the  software  to  overcome  another  serious  problem,  that  of  backup. 
Orailoglu  [20],  in  his  masters  thesis,  created  an  experimental 
implementation  of  a  multirate  method  using  a  separate  integrator  for  each 
equation.   Each  integrator  was  actually  a  form  of  DIFSUB.   An  "event 
driven"  approach  was  used.   The  basic  idea  is  to  examine  the  estimated 
stepsizes  for  all  equations,  and  to  choose  that  equation  (and  corresponding 
variable)  for  integration  which  will  progress  to  the  earliest  point  in 
time.   The  variable  is  then  integrated,  and  the  estimated  stepsize  for  the 
next  step  in  that  variable  calculated.   Each  mesh  point  in  each  equation 
can  be  thought  of  as  an  "event"  from  a  simulation  viewpoint.   The 
simulation  proceeds  by  advancing  to  the  next  event  to  occur.   The  problem 
with  this  approach  is  that  the  stepsizes  for  subsequent  steps  can  only  be 
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estimated.   If  when  the  step  is  actually  taken  it  turns  out  to  be  too 
large,  it  must  be  reduced.   The  difficulty  arises  when  one  variable  is 
being  integrated  with  a  very  small  stepsize  and  another  with  a  very  large 
one  that  turns  out  too  be  much  too  large.   For  example,  suppose  that  y  and 
z  are  being  integrated  with  stepsizes  h  and  lOOh  repectively.   If  they  are 
both  currently  known  at  time  0,   99  or  100  integrations  of  y  will  be 
performed  before  one  integration  of  z  is  attempted.  When  we  attempt  to 
integrate  z  and  find  that  the  estimated  stepsize  is  much  too  large,  we  will 
reduce  the  stepsize  and  find  ourselves  trying  to  integrate  z  to  a  much 
earlier  point  in  time  than  the  current  t  value  of  y.   Consequently,  we  will 
not  have  any  values  of  y  on  which  to  base  an  interpolation  for  the 
evaluation  of  the  derivative  of  z.   This  problem  may  not  be  as  serious  as 
the  fact  that  y  has  been  integrated  using  extrapolated  values  of  z  over  an 
interval  that  is  too  large  for  an  integration  step.   Since  the  error  bound 
for  interpolation  or  extrapolation  has  the  same  general  form  as  the  error 
bound  for  integration,  there  is  a  danger  that  the  integration  of  y  will 
have  been  contaminated  with  interpolation  errors  from  z.   Orailoglu  tried 
various  techniques  to  avoid  this  problem.   One  method  was  to  use  "pre- 
integration"  where  every  step  was  done  twice,  with  the  first  step  used  to 
get  a  better  estimate  of  the  stepsize.   Even  though  this  was  expensive,  the 
total  number  of  function  evaluations  was  still  less  than  in  a  conventional 
method.   (Here,  a  function  evaluation  is  an  evaluation  of  a  single  f  .) 
Another  mechanism  was  tried  with  some  success .   We  monitored  the  coupling 
of  errors  from  variable  to  variable  by  examining  the  Jacobian  of  the 
system.   Unfortunately,  this  approach  required  knowledge  of  the  Jacobian 
even  for  non-stiff  systems.   We  were  able  to  "predict"  when  an  estimated 
stepsize  had  to  be  reduced  by  this  technique,  but  the  cost  was  high.   We 
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have  started  to  investigate  a  different  approach  based  on  the  observation 
that  the  coupling  from  a  rapidly  varying  component  to  a  slowly  varying 
component  must  be  weak  (or  the  latter  would  not  be  slowly  varying)  . 
Therefore,  it  may  be  possible  to  to  integrate  the  slowly  changing 
components  over  their  large  stepsizes  first,  their  values  at  the  end  of  the 
interval  being  relatively  independent  of  the  fast  components.   Then,  the 
fast  components  can  be  integrated,  using  more  accurate  interpolation  for 
the  intermediate  values  of  the  slow  components  rather  than  extrapolation. 
The  process  looks  like  a  top-down,  successive  refinement  process  because  it 
can  be  implemented  in  the  following  way: 

(i)    Initially,  assume  all  variables  are  known  at  the  same  time  value, 
say  t  =  0. 

(ii)    Select  the  set  of  variables  known  at  the  earliest  time.   (Initially, 
this  is  the  complete  set.)  Integrate  these  equations  with  the 
largest  of  the  estimated  stepsizes  of  any  of  the  variables. 

(iii)   Check  the  error  of  each  variable  individually. 

(iv)   Reject  the  steps  on  the  variables  with  large  errors,  and  halve  their 
suggested  stepsizes. 

(v)    Repeat  the  whole  process  from  step  (ii) • 

This  technique  performs  a  compound  step  whose  length  is  the  largest 
stepsize  for  which  one  equation  can  be  integrated  successfully.   It  does 
successive  halving  of  other  stepsizes,  and  reintegrates  them  until  the 
stepsize  is  small  enough.   Although  this  may  lead  to  many  wasted  steps,  the 
fact  that  the  binary  subdivision  reaches  the  stepsize  in  a  logarithmic 
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number  of  subdivisions  means  that  the  number  of  failed  steps  never  exceeds 
the  number  of  successful  steps.   While  50%  efficiency  may  seem  poor,  the 
name  of  the  game  is  to  reduce  the  number  of  integrations  in  the  slow 
components  to  such  an  extent  that  we  are  ahead.   When  applied  to  equations 
arising  from  the  method  of  lines,  this  approach  automatically  refines  the 
mesh  in  regions  of  rapid  transition. 

There  are,  however,  many  unanswered  questions  with  this  technique. 
For  example,  the  choice  of  the  outermost  large  step  is  critical,  so  a 
theory  to  guide  its  selection  is  needed.   Error  estimates  are  also 
critical.   If  an  estimate  is  small  for  a  particular  variable,  its  step  will 
not  be  subdivided  further,  even  though  subdivision  of  steps  for  other 
variables  may  cause  so  much  change  that  the  integration  of  the  first 
variable  was  invalid. 

^.^_.   Low  accuracy  methods 

Current  ODE  methods  are  based  on  the  concept  of  polynomial  order. 
This  matches  the  numerical  solution  to  the  true  solution  in  an  asymptotic 
sense.   Unfortunately,  for  large  stepsizes  several  problems  arise.   The 
first  is  that  error  estimates  are  unreliable  because  they  also  are  based  on 
asymptotic  theory.   The  second  is  that  the  actual  methods  are  inefficient 
and  display  poor  stability  properties.   Low  accuracy  methods  are  needed  in 
a  number  of  applications,  in  particular,  real-time  simulation  [15]  and  the 
time  integration  of  PDEs.   In  real-time  integration,  one  of  the  most 
important  properties  is  the  stability  correspondence  between  the  actual 
system  and  the  numerical  model  because  the  simulation  is  frequently  being 
done  to  test  the  "handling  properties"  of  a  proposed  device.   The  method 
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should  be  stable  or  unstable  as  the  actual  device  Is.   The  same 
characteristics  are  probably  desirable  for  applications  to  PDEs .   Since  we 
are  talking  about  a  large  stepsize,  polynomial  accuracy  is  probably  not 
important.   Some  form  of  minimax  approximation  is  more  appropriate  (as  it 
is  in  approximation  over  an  interval  rather  than  in  the  neighborhood  of  a 
point).   For  example,  if  we  consider  the  class  of  differential  equations 

y'  =  xy 

where  X   is  known  to  lie  in  some  region  D  of  the  complex  plane,  we  can  find 
the  Runge-Kutta  or  multistep  method  with  a  given  number  of  stages  or  steps 
that  gives  the  best  approximation  to  the  class  of  equations  in  the  minimax 
sense.   The  approximation  also  depends  on  the  stepsize  h,  so  the  actual 
maximum  error  depends  on  the  stepsize.   In  fact,  it  can  be  seen  that  this 
dependency  is  monotonic,  so  the  step  can  be  chosen  when  the  region  D  is 
given  and  the  error  tolerance  is  specified.   The  general  approach  would  be 
to  assume  a  form  for  the  class  of  problems  with  some  parameters.   Formally, 
we  would  have  a  class  of  problems  C(p)  where  the  parameters  p  are  assumed 
to  lie  in  a  region  D  of  the  parameter  space.   A  class  of  methods  M(h,a) 
would  be  considered,  where  a  is  a  set  of  method  coefficients,  and  h  is  the 
stepsize.   For  a  given  stepsize  h,  the  coefficients  should  be  picked  to 
minimize  the  worst  case  error  of  the  class  of  problems.   Thus,  if  the  error 
for  problem  C(p)  using  method  M(h,a)  is  E(p,h,a),  we  would  pick  a  for  a 
given  h  to  minimize  the  worst  case  error  to  get 


E(h)  =  min  max  E(p,h,a) 
a  p  in  D 

The  stepsize  would  be  picked  so  that  E(h)  is  sufficiently  small,  that  is, 
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so  that  E(h)  <  e.   An  automatic  program  might  function  by  estimating  the 
regions  in  which  the  parameters  lie,  and  using  that  information  to  reduce 
the  region  over  which  the  minimax  approximation  is  taken.   The  utility  of 
an  approach  of  this  form  is  dependent  on  determining  some  simple 
approximations  to  the  minimax  solutions.   Clearly  one  cannot  hope  to  solve 
a  complex  minimization  problem  at  each  integration  step! 

4^.   Conclusion 

This  paper  has  surveyed  a  number  of  recent  advances,  many  of  which  are 
very  important,  either  offering  major  advances  in  the  theory  or  suggesting 
ways  of  handling  currently  very  difficult  problems.   It  has  also  examined 
some  problems  that  remain  very  challenging,  and  these  are  problems  whose 
solution  can  have  a  significant  impact  on  the  solution  of  PDE  problems. 
Some  of  these  problems  have  appeared  very  recently,  which  suggests  that 
there  will  continue  to  be  new  problems  to  solve  in  this  field  for  some 
time,  and  that  continued  effort  can  reap  worthwhile  rewards. 
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